1. Introduzione

Il presente progetto si propone di applicare diverse tecniche di clustering al dataset banknote al fine di identificare gruppi omogenei presenti nei dati e confrontare l’efficacia dei vari metodi. L’analisi intende mettere in pratica le metodologie apprese nel modulo di Data Analysis & Statistical Learning, in particolare, verranno trattati i seguenti algoritmi di clustering: Clustering Gerarchico Agglomerativo e Divisivo, Clustering Partizionale (K-means e K-medoids) e Clustering Model-based (Mixture of Gaussian).

Tali algoritmo sono stati valutati su tre diverse tipologie di dataset: - Dataset Completo: utilizzo di tutte le variabili disponibili. - Dataset Ridotto: uso esclusivo delle variabili “Bottom” e “Diagonal”, individuate come le più discriminative. - Dataset Trasformato: applicazione della PCA per ridurre le variabili a componenti che sintetizzano l’informazione.

Per ciascuna configurazione, sono state valutate le performance dei metodi di clustering mediante metriche interne ed esterne. Il documento è organizzato nei seguenti capitoli:

  1. Introduzione: Presentazione del progetto, degli obiettivi e del contesto.
  2. Descrizione del Dataset e Analisi Esplorativa: Esplorazione preliminare dei dati, analisi descrittiva e operazioni di preprocessing, tra cui la costruzione del dataset trasformato mediante la tecnica PCA.
  3. Applicazione delle Tecniche di Clustering: Implementazione dei metodi di clustering (gerarchico, partizionale e model-based).
  4. Validazione e Visualizzazione dei Risultati: Confronto dei risultati mediante metriche di validazione e visualizzazione dei cluster (con l’ausilio della PCA).
  5. Conclusioni e Sviluppi Futuri: Sintesi dei risultati ottenuti e prospettive per ulteriori analisi.

2. Descrizione e analisi del dataset

In questa sezione esamineremo il dataset banknote per comprenderne la struttura, le caratteristiche e le relazioni tra le variabili. Il dataset proviene dal pacchetto mclust e contiene misure che possono essere utilizzate per distinguere tra banconote autentiche e false (o per altri scopi diagnostici). Il dataset contiene le seguenti variabili:

Di seguito carichiamo tutte le librerie necessarie per il progetto.

library("mclust")
library("reshape2")
library("ggplot2")
library("GGally")
library("dplyr")
library("clustertend")
library("FNN")
library("FactoMineR")
library("factoextra")
library("caret")
library("pheatmap")
library("clusterCrit")
library("NbClust")
library("fpc")
library("scatterplot3d")
library("corrplot")
library("cluster")
library("pROC")
library("gridExtra")

Carichiamo il dataset.

data(banknote)

Visualizziamo le prime righe del dataset e la sua struttura.

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4
str(banknote)
## 'data.frame':    200 obs. of  7 variables:
##  $ Status  : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Length  : num  215 215 215 215 215 ...
##  $ Left    : num  131 130 130 130 130 ...
##  $ Right   : num  131 130 130 130 130 ...
##  $ Bottom  : num  9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
##  $ Top     : num  9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 11 10 ...
##  $ Diagonal: num  141 142 142 142 142 ...

Il dataset è composto da 200 osservazioni e 7 variabili, tutte di tipo numerico ad eccezione di Status che è una variabile binaria che può assumere uno tra i seguenti valori: counterfeit e genuine.

Visualizziamo un riepilogo statistico delle variabili.

summary(banknote)
##          Status        Length           Left           Right      
##  counterfeit:100   Min.   :213.8   Min.   :129.0   Min.   :129.0  
##  genuine    :100   1st Qu.:214.6   1st Qu.:129.9   1st Qu.:129.7  
##                    Median :214.9   Median :130.2   Median :130.0  
##                    Mean   :214.9   Mean   :130.1   Mean   :130.0  
##                    3rd Qu.:215.1   3rd Qu.:130.4   3rd Qu.:130.2  
##                    Max.   :216.3   Max.   :131.0   Max.   :131.1  
##      Bottom            Top           Diagonal    
##  Min.   : 7.200   Min.   : 7.70   Min.   :137.8  
##  1st Qu.: 8.200   1st Qu.:10.10   1st Qu.:139.5  
##  Median : 9.100   Median :10.60   Median :140.4  
##  Mean   : 9.418   Mean   :10.65   Mean   :140.5  
##  3rd Qu.:10.600   3rd Qu.:11.20   3rd Qu.:141.5  
##  Max.   :12.700   Max.   :12.30   Max.   :142.4

Notiamo che la colonna Status è perfettamente bilanciata.

Controlliamo se ci sono valori mancanti nelle colonne.

colSums(is.na(banknote))
##   Status   Length     Left    Right   Bottom      Top Diagonal 
##        0        0        0        0        0        0        0

Come possiamo vedere nessuna colonna presenta valori mancanti. Per avere una prima idea delle relazioni tra le variabili, creiamo una matrice di scatterplot. Utilizziamo la colonna Status per colorare i punti.

ggpairs(banknote[,-1],upper = list(continuous = "density", combo = "box_no_facet"),
        lower = list(continuous = "points", combo = "dot_no_facet"),aes(colour=banknote$Status))

TODO: sintetizzare Lungo la diagonale sono visibili le distribuzioni univariate delle variabili. Si nota una sovrapposizione parziale tra i due gruppi di colore (Status), ma anche zone in cui uno dei due colori prevale. Nei grafici di densità bivariata (metà diagonale superiore della matrice) e negli scatter plot (metà diagonale inferiore della matrice), si osservano alcune regioni in cui i punti di uno stesso colore tendono a concentrarsi, suggerendo cluster parzialmente distinti. Nonostante la sovrapposizione in diverse aree, si intravedono regioni in cui un gruppo prevale. Ciò suggerisce che un algoritmo di classificazione o di clustering potrebbe distinguere parzialmente le due classi, anche se non in modo perfetto con una singola coppia di variabili. È importante sottolineare che questo grafico mostra solo le relazioni a coppie. Ciò significa che l’analisi si limita a considerare due variabili alla volta, mentre una clusterizzazione netta potrebbe emergere solo quando si considerano contemporaneamente più variabili. In altre parole, anche se alcune coppie non evidenziano una separazione chiara, combinando più dimensioni si potrebbe scoprire una struttura clusterizzata più marcata. Il fatto che in alcune proiezioni (coppie di variabili) si vedano zone di separazione o contorni distinti suggerisce che esistono sottostrutture nei dati. Notiamo infine che le variabili Bottom e Diagonal sembrerebbero separare in modo migliore il dataset rispetto alle altre variabili.

D’ora in avanti, per ciascun algoritmo di clustering che proveremo, verranno effettuate tre diverse analisi:

  1. Dataset completo: utilizzando tutte le variabili disponibili nel dataset.

  2. Dataset ridotto: considerando esclusivamente le variabili Diagonal e Bottom.

  3. Dataset trasformato: Applicando l’Analisi delle Componenti Principali (PCA) per ridurre la dimensionalità.

Queste tre approcci permetteranno di valutare l’impatto delle diverse selezioni di variabili e della riduzione dimensionale sulle prestazioni e sui risultati degli algoritmi di clustering.

2.1. Valutazione della tendenza al clustering

Prima di applicare gli algoritmi di clustering, è fondamentale verificare se il dataset presenta una naturale propensione alla formazione di gruppi. Utilizziamo specifiche misure statistiche per determinare se i dati sono strutturati in cluster o se, al contrario, sono distribuiti in modo casuale. Questi metodi permettono di valutare se il clustering è effettivamente significativo e utile per l’analisi del dataset.

2.1.1. Hopkins statistics

L’Hopkins Statistic è una misura utilizzata per valutare la tendenza di un dataset alla clusterizzazione, confrontando la sua struttura con quella di un insieme di punti generati casualmente. Assume valori compresi tra 0 e 1, più il valore è vicino allo 0 più il dataset è strutturato in cluster.

banknote_num <- banknote[, -1]

banknote_num <- scale(banknote_num)

set.seed(987)

hopkins_stat <- hopkins(banknote_num, n = nrow(banknote_num) - 1)

random_df <- apply(banknote_num, 2, function(x){runif(length(x), min(x), max(x))})
random_df <- as.data.frame(random_df)

random_df <- scale(random_df)

set.seed(987)

hopkins_stat_rand <- hopkins(random_df, n = nrow(random_df)-1)

cat("Hopkins Statistics on Banknote dataset:", hopkins_stat$H, "\n")
## Hopkins Statistics on Banknote dataset: 0.2642237
cat("Hopkins Statistics on Random dataset:", hopkins_stat_rand$H, "\n")
## Hopkins Statistics on Random dataset: 0.4926061

Un valore di Hopkins pari a 0.2519264, essendo inferiore a 0.5, indica una tendenza alla clusterizzazione. In altre parole, le distanze tra i punti reali sono molto più ridotte rispetto a quelle dei punti generati casualmente, suggerendo la presenza di gruppi ben definiti nel dataset. Questo risultato supporta l’ipotesi che il dataset sia strutturato in cluster distinti.

2.1.2. Visual Assessment of cluster Tendency (VAT)

L’algoritmo VAT produce una matrice di dissimilarità ordinata, ottenuta riorganizzando le righe e le colonne della matrice delle distanze per evidenziare eventuali blocchi di punti “ravvicinati” (ossia potenziali cluster). Visualizziamo la matrice con l’intento di verificare se si formano lungo la diagonale secondaria del grafico, dei blocchi rossi che rappresentano regioni di osservazioni che condividono una forte similarità.

p1 <- fviz_dist(dist(banknote), show_labels = FALSE) + labs(title = "Banknote data")
p2 <- fviz_dist(dist(random_df), show_labels = FALSE) + labs(title = "Random data")

grid.arrange(p1, p2, ncol = 2, widths = c(6, 6))

Nei grafici VAT il rosso indica alta similarità (ossia bassa dissimilarità), mentre il blu indica bassa similarità (ossia alta dissimilarità). Nel grafico del dataset banknote, notiamo che lungo la diagonale sono presenti regioni rosse che indicano dati simili tra loro, suggerendo la presenza di cluster. Al contrario, il dataset randomico appare privo di struttura.

Possiamo concludere che il dataset Banknote ha una tendezza alla clusterizzazione, procediamo perciò con il calcolo della PCA e successivamente alla clusterizzazione.

2.2. Riduzione della dimensionalità con PCA

In questa sezione applichiamo la PCA sul dataset, considerando solo le colonne numeriche. Anche se le variabili sono tutte misurate in millimetri e quindi nella stessa scala, applichiamo comunque la standardizzazione per evitare che differenze nella dispersione (variabili con varianze maggiori) influenzino in modo sproporzionato i risultati della PCA.

banknote_num <- banknote[, -1]
head(banknote_num)
##   Length  Left Right Bottom  Top Diagonal
## 1  214.8 131.0 131.1    9.0  9.7    141.0
## 2  214.6 129.7 129.7    8.1  9.5    141.7
## 3  214.8 129.7 129.7    8.7  9.6    142.2
## 4  214.8 129.7 129.6    7.5 10.4    142.0
## 5  215.0 129.6 129.7   10.4  7.7    141.8
## 6  215.7 130.8 130.5    9.0 10.1    141.4
pca_result <- PCA(banknote_num, scale = TRUE, graph = FALSE)

Estraiamo gli autovalori e calcoliamo sia la varianza spiegata per ciascuna componente sia la varianza cumulativa. Queste informazioni ci serviranno per decidere il numero ottimale di componenti principali da mantenere.

eigenvalues <- get_eigenvalue(pca_result)[,1]
explained_variance <- eigenvalues / sum(eigenvalues)
cumulative_variance <- cumsum(explained_variance)

print(eigenvalues)
##     Dim.1     Dim.2     Dim.3     Dim.4     Dim.5     Dim.6 
## 2.9455582 1.2780838 0.8690326 0.4497687 0.2686769 0.1888799
print(explained_variance)
##      Dim.1      Dim.2      Dim.3      Dim.4      Dim.5      Dim.6 
## 0.49092637 0.21301396 0.14483876 0.07496145 0.04477948 0.03147998
print(cumulative_variance)
##     Dim.1     Dim.2     Dim.3     Dim.4     Dim.5     Dim.6 
## 0.4909264 0.7039403 0.8487791 0.9237405 0.9685200 1.0000000

Abbiamo ora a disposizione gli autovalori e le informazioni sulla varianza, utili per valutare l’importanza di ciascuna componente principale.

Utilizziamo la Kaiser Rule, una tecnica utilizzata per determinare il numero di componenti principali da mantenere quando si esegue l’analisi delle componenti principali. La regola si basa sull’idea di conservare solo le componenti con una varianza spiegata maggiore di 1. Questo criterio si applica al valore proprio (eigenvalue) delle componenti principali.

kaiser_components <- sum(eigenvalues > 1)
cat("Numero di componenti selezionati con la Kaiser rule:", kaiser_components, "\n")
## Numero di componenti selezionati con la Kaiser rule: 2

Otteniamo 2 come numero di componenti principali da mantenere, ciò significa che le due componenti principali con valori propri superiori a 1 spiegano una quantità di varianza maggiore di quella che verrebbe spiegata da una singola variabile originale del dataset.

Un ulteriore criterio è quello di selezionare il numero minimo di componenti la cui somma di varianza spiegata (varianza comulativa), raggiunge almeno l’80% del totale. In questo modo garantiamo una buona rappresentazione dei dati.

variance_threshold <- min(which(cumulative_variance >= 0.80))

cat("Numero di componenti principali che spiegano l'80% della varianza comulativa:", variance_threshold, "\n")
## Numero di componenti principali che spiegano l'80% della varianza comulativa: 3
fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 50))

Nel grafico vediamo la varianza spiegata da ciascuna componente. Con le prime tre componenti principali riusciamo a spiegare l’84.9% della varianza cumulativa.

Vediamo infine un ultimo metodo per scegliere il numero di componenti principali: l’elbow method. Mediante un grafico scree plot visualizziamo gli autovalori per ciascuna componente e individuiamo “il gomito”, ovvero il punto nel grafico in cui la diminuzione degli autovalori rallenta. Tale punto può essere interpretato come un’indicazione del numero ottimale di componenti. Includiamo nel grafico anche due linee verticali per evidenziare i criteri della varianza cumulativa (rosso) e della Kaiser Rule (blu).

scree_plot <- data.frame(PC = 1:length(eigenvalues), Eigenvalue = eigenvalues)

ggplot(scree_plot, aes(x = PC, y = Eigenvalue)) +
    geom_point() +
    geom_line() +
    geom_vline(xintercept = variance_threshold, linetype = "dashed", color = "red") +
    geom_vline(xintercept = kaiser_components, linetype = "dotted", color = "blue") +
    labs(title = "Scree Plot for PCA", x = "Principal Component", y = "Eigenvalue") +
    theme_minimal()

Nel grafico il gomito sembra collocarsi intorno alla quarta componente. I tre criteri danno perciò risultati differenti, per mantenere un equilibrio tra interpretabilità e preservazione della varianza, si potrebbe optare per 3 componenti.
Di seguito viene costruito il nuovo dataset in cui le variabili sono sostituite con le prime 3 PCA, che verrà usato più avanti nelle analisi.

pca_data <- as.data.frame(pca_result$ind$coord[, 1:3])
colnames(pca_data) <- c("PCA1", "PCA2", "PCA3")
banknote_pca <- cbind(pca_data, Status = banknote$Status)
banknote_pca <- banknote_pca %>% dplyr::select(Status, everything())
head(banknote_pca)
##    Status       PCA1        PCA2       PCA3
## 1 genuine  1.7474012  1.65082829  1.4237611
## 2 genuine -2.2743177 -0.53879328  0.5326484
## 3 genuine -2.2774015 -0.10767707  0.7174149
## 4 genuine -2.2835546 -0.08765431 -0.6056336
## 5 genuine -2.6321283  0.03919590  3.1963847
## 6 genuine  0.7584073  3.08874513  0.7864804

Per comprendere meglio il contributo delle variabili alle componenti principali, realizziamo una mappa di correlazione. Questo ci aiuta a capire come le variabili si relazionano alle componenti principali.

var <- get_pca_var(pca_result)

fviz_pca_var(pca_result, col.var = "contrib",
            gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Questo passaggio ci offre una visione immediata di come ogni variabile contribuisce alle componenti principali e della loro interrelazione. Osserviamo che le variaibili che appaiono con colori più accesi sono Length che contribuisce molto a Dim.2, e le variabili Diagonal, Left e Right con un contributo importante a Dim.1 seppure opposto.

Approfondiamo l’analisi mostrando i contributi delle variabili attraverso diversi metodi grafici, in modo da identificare quali variabili influenzano maggiormente ciascuna componente.

corrplot(var$contrib, is.corr=FALSE)

L’analisi dei contributi evidenzia il peso relativo di ciascuna variabile nelle componenti principali, utile per interpretare i risultati della PCA. Si osserva che Length e Top hanno contributi particolarmente elevati rispettivamente a Dim.2 e Dim.3.

Approfondiamo come hanno contribuito le variabili per la prima componente principale.

fviz_contrib(pca_result, choice = "var", axes = 1)

Come era visibile già dal grafico precedente, le variabili che hanno contribuito maggiormente alla prima componente principale sono Diagonal, Right, Left.

Similmente al passaggio precedente, visualizziamo i contributi delle variabili per la seconda componente principale.

fviz_contrib(pca_result, choice = "var", axes = 2)

Notiamo che la variabile Length ha contribuito più del 60% alla seconda componente principale.

3. Clustering gerarchico

3.1 Clustering Gerarchico Agglomerativo

Il Clustering Gerarchico Agglomerativo (HAC) è una tecnica di clustering che costruisce una gerarchia di cluster attraverso un approccio “bottom-up”. Inizia considerando ogni punto dati come un cluster separato e, successivamente, unisce iterativamente i cluster più simili fino a formare un unico cluster che racchiude tutti i dati. Questo metodo non richiede la specifica del numero di cluster a priori e produce una rappresentazione ad albero, nota come dendrogramma, che facilita la visualizzazione delle relazioni tra i cluster.

Per applicare il clustering gerarchico agglomerativo al nostro dataset, seguiamo questi passaggi:

  1. Scaling: Standardizziamo i dati.

  2. Calcolo della Matrice di Distanza: Determiniamo le distanze tra tutte le coppie di punti nel dataset. Una misura comune è la distanza euclidea.

  3. Applicazione del clustering agglomerativo: Utilizziamo la funzione hclust() per eseguire il clustering gerarchico e generare il dendrogramma.

  4. Taglio del Dendrogramma: Decidiamo il numero ottimale di cluster utilizzando varie tecniche, e tagliamo il dendrogramma al livello appropriato.

  5. Valutazione: Valutiamo i risultati del clustering utilizzando varie metriche.

3.1.1. Dataset completo

Scaliamo il dataset e calcoliamo la relativa matrice di dissimilatità mediante distanza euclidea.

df <- scale(banknote_num)
distanze <- dist(df, method = "euclidean")
diss_matrix <- as.matrix(distanze)

Applichiamo il clustering con metodo di linkage ward che generalmente insieme all’average linkage, performa meglio.

hclust_clustering <- hclust(distanze, method = "ward.D2")

Il dendrogramma mostra come le osservazioni si uniscono in cluster, con l’altezza dei rami che riflette la loro dissimilarità. Analizzando il dendrogramma, possiamo decidere il numero di cluster, tagliando l’albero a un’altezza che separa i gruppi in modo significativo.

3.1.1.1. Scelta del miglior numero di cluster k

fviz_dend(hclust_clustering, show_labels = FALSE, as.ggplot = TRUE)

Notiamo un ampio salto a circa 30 suggerisce una buona separazione in 2 cluster, mentre un altro salto meno marcato (intorno a 15) ne suggerisce 3. Per confermare il numero ottimale di cluster, si prosegue con metodi alternativi.

Osserviamo l’andamento del Total Within Sum of Squares (WSS) rispetto a k, cercando di individuare il punto in cui la curva inizia a “piegarsi” più marcatamente, cioè dove l’aggiunta di ulteriori cluster non porta a una riduzione significativa del WSS. Di seguito il grafico Elbow.

fviz_nbclust(df, FUN = hcut, method = "wss") +
    geom_vline(xintercept = 2, linetype = 2) +
    labs(subtitle = "Elbow method")

Tra k = 1 e k = 2, la diminuzione del WSS è molto forte, tra k = 2 e k = 3 c’è ancora un calo del WSS ma meno netto. Oltre k = 3, la pendenza inizia a ridursi gradualmente. Dal grafico, sembra che il calo principale avvenga per k = 2, mentre oltre si riduce più lentamente.

Analizziamo ora il grafico dell’Average Silhouette Width.

fviz_nbclust(df, FUN = hcut, method = "silhouette") +
    labs(subtitle = "Silhouette method")

L’indice di silhouette risulta più alto in corrispondenza di k = 2. Questo suggerisce che la separazione tra i cluster e la coesione interna sia massima quando si dividono i dati in due gruppi. Da k = 3 l’indice diminuisce senza tornare ai valori di picco.

Osserviamo il grafico del Gap Statistic che confronta la dispersione all’interno dei cluster ottenuta dal clustering, con quella attesa in un dataset generato casualmente (quindi privo di gruppi).

set.seed(123)
fviz_nbclust(df, FUN = hcut, method = "gap_stat", nboot = 500, verbose = FALSE) +
    labs(subtitle = "Gap statistic method")

Il grafico del Gap Statistic indica che quando il dataset viene diviso in 3 gruppi la struttura di clustering risulta migliore rispetto a quella attesa in un dataset casuale.

Utilizziamo infine il pacchetto NbClust che valuta il numero ottimale di cluster utilizzando numerosi indici interni.

pdf(NULL)
nb <- NbClust(df, distance = "euclidean", min.nc = 2, max.nc = 10, method = "ward.D2")
## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 
## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 7 proposed 2 as the best number of clusters 
## * 10 proposed 3 as the best number of clusters 
## * 1 proposed 4 as the best number of clusters 
## * 1 proposed 7 as the best number of clusters 
## * 1 proposed 8 as the best number of clusters 
## * 3 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************
dev.off()
## png 
##   2

La decisione finale basata sulla majority rule suggerisce 3 cluster, indicando che, complessivamente, questa soluzione offre una migliore separazione e coesione interna rispetto alle altre soluzioni.

Complessivamente la metà dei metodi analizzati (gap statistic, nbclust) converge su una soluzione a 3 cluster, mentre l’altra metà (elbow, e silhouette) converge a 2 cluster. Pertanto non possiamo optare per una scelta di “maggioranza”, decidiamo di scegliere un numero di cluster pari a 2 poichè il salto più marcato nel dendrogramma indica una forte separazione naturale in 2 gruppi.

Procediamo con il taglio del dendogramma in modo da ottenere 2 cluster. Successivamente visualizziamo una tabella di frequenza che mostra il numero di osservazioni assegnate a ciascun cluster, permettendoci di verificare la distribuzione dei dati tra i 3 cluster ottenuti.

cluster_cut <- cutree(hclust_clustering, k = 2)

table(cluster_cut)
## cluster_cut
##   1   2 
## 102  98
fviz_dend(hclust_clustering, k = 2,
        cex = 0.5,
        k_colors = c("#2E9FDF", "#E7B800"),
        color_labels_by_k = TRUE,
        rect = TRUE
)

3.1.1.2. Validazione

Valutiamo ora la qualità dei risultati del nostro algoritmo di clustering mediante quattro metriche: Average Silhouette Width, Dunn Index, Corrected Rand Index e Meila’s Variation of Information Index.

Iniziamo con la metrica della silhouette media (calcolata come la media di s(i) su tutte le osservazioni), che fornisce un’indicazione complessiva della qualità del clustering: si auspica un valore il più vicino ad 1 possibile.

sil <- silhouette(cluster_cut, distanze)
fviz_silhouette(sil, palette = "jco", ggtheme = theme_classic())
##   cluster size ave.sil.width
## 1       1  102          0.39
## 2       2   98          0.36

Un valore di 0.37 suggerisce che la struttura dei cluster non è particolarmente forte: c’è una separazione moderata, ma c’è anche una certa sovrapposizione tra i cluster.

Analizziamo ora il Dunn Index, che valuta la separazione tra i cluster rispetto alla loro compattezza, cioè quanto i cluster siano ben separati e internamente compatti. Se i dati contengono cluster compatti e ben separati, ci si aspetta che il diametro dei cluster sia piccolo e che la distanza tra i cluster sia grande. Pertanto l’indice di Dunn deve essere massimizzato.

dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn

print(paste("Dunn Index:", dunn_index))
## [1] "Dunn Index: 0.210791490416207"

Un Dunn Index di 0.210 indica che i cluster hanno una separazione piuttosto bassa rispetto alla loro coesione interna, suggerendo possibili problemi di sovrapposizione.

Analizziamo il Corrected Rand Index che misura quanto bene il clustering ottenuto corrisponde a una partizione di riferimento. Nel nostro caso utilizziamo la colonna Status del dataset. Ed infine il Meila’s Variation of Information (VI) Index che misura quanto due cluster differiscono tra loro in termini di entropia. A differenza del Corrected Rand Index, che valuta la similarità, il VI Index misura la distanza tra due partizioni (va perciò minimizzato).

status <- as.numeric(banknote$Status)
clust_stats <- cluster.stats(d = dist(df), status, cluster_cut)
print(paste("Corrected Rand Index:", clust_stats$corrected.rand))
## [1] "Corrected Rand Index: 0.960200080403878"
print(paste("Meila's Index:", clust_stats$vi))
## [1] "Meila's Index: 0.0982391266151992"

Un Corrected Rand Index di 0.9602 è molto elevato, il che indica che il clustering ottenuto corrisponde molto bene alla partizione di riferimento. Inoltre, un Information Index di 0.098 è basso, suggerendo che la differenza informativa tra le partizioni è minima.

Concludendo ci potrebbe essere sovrapposizione tra i cluster, ma senza compromettere eccessivamente la qualità dell’assegnazione.

3.2.1.3. Ottimizzazione

Ottimizziamo il clustering esplorando diverse combinazioni di metodi di linkage e numero di cluster. Valutiamo ciascuna combinazione mediante le metriche viste prima. Successivamente, assegniamo un ranking per ciascuna metrica (invertendo il criterio per quelle da massimizzare) e sommiamo i punteggi. Infine, ordiniamo i risultati e selezioniamo la combinazione con il punteggio complessivo più basso come miglior clustering.

evaluate_clustering <- function(df, distanze, status, linkage_methods = c("ward.D2", "single", "complete", "average", "centroid"), min_k = 2, max_k = 10) {
    results_all <- data.frame(Method = character(),
                            k = numeric(),
                            Silhouette = numeric(),
                            Dunn = numeric(),
                            Corrected_Rand = numeric(),
                            VI = numeric(),
                            stringsAsFactors = FALSE)

    for(method in linkage_methods) {
        hclust_obj <- hclust(distanze, method = method)
        
        for(k in min_k:max_k) {
            cluster_cut <- cutree(hclust_obj, k = k)
            sil <- silhouette(cluster_cut, distanze)
            avg_sil <- mean(sil[, 3])
            dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
            clust_stats <- cluster.stats(distanze, status, cluster_cut)
            corrected_rand <- clust_stats$corrected.rand
            vi_val <- clust_stats$vi
            
            results_all <- rbind(results_all, data.frame(Method = method,
                                                        k = k,
                                                        Silhouette = avg_sil,
                                                        Dunn = dunn_index,
                                                        Corrected_Rand = corrected_rand,
                                                        VI = vi_val,
                                                        stringsAsFactors = FALSE))
        }
    }
    return(results_all)
}

results_all <- evaluate_clustering(df, distanze, status)

best_silhouette <- results_all %>% group_by(Method) %>% filter(Silhouette == max(Silhouette))
print(best_silhouette)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete     2      0.346 0.125          0.591 0.566 
## 4 average      2      0.333 0.341          0     0.718 
## 5 centroid     2      0.333 0.341          0     0.718
best_dunn <- results_all %>% group_by(Method) %>% filter(Dunn == max(Dunn))
print(best_dunn)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete    10      0.186 0.179          0.268 1.45  
## 4 average      2      0.333 0.341          0     0.718 
## 5 centroid     2      0.333 0.341          0     0.718
best_corrected_rand <- results_all %>% group_by(Method) %>% filter(Corrected_Rand == max(Corrected_Rand))
print(best_corrected_rand)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2     0.375  0.211       0.960    0.0982
## 2 single       8    -0.0124 0.225       0.000320 0.896 
## 3 complete     3     0.327  0.138       0.688    0.596 
## 4 average      5     0.316  0.222       0.951    0.140 
## 5 centroid     3     0.228  0.298       0.000101 0.742
best_vi <- results_all %>% group_by(Method) %>% filter(VI == min(VI))
print(best_vi)
## # A tibble: 5 × 6
## # Groups:   Method [5]
##   Method       k Silhouette  Dunn Corrected_Rand     VI
##   <chr>    <int>      <dbl> <dbl>          <dbl>  <dbl>
## 1 ward.D2      2      0.375 0.211          0.960 0.0982
## 2 single       2      0.333 0.341          0     0.718 
## 3 complete     2      0.346 0.125          0.591 0.566 
## 4 average      5      0.316 0.222          0.951 0.140 
## 5 centroid     2      0.333 0.341          0     0.718
results_ranked <- results_all %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall <- results_ranked %>% arrange(total_rank) %>% head(1)
print(best_overall)
##    Method k Silhouette      Dunn Corrected_Rand         VI rank_sil rank_dunn
## 1 ward.D2 2  0.3748428 0.2107915      0.9602001 0.09823913        1        21
##   rank_corrected_rand rank_vi total_rank
## 1                   1       1         24

Il risultato mostra che il clustering ottimale secondo queste metriche è Ward.D2 con 2 cluster, che è la combinazione che abbiamo testato precedentemente. La combinazione scelta restituisce i valori migliori possibili per le metriche Silhouette, Corrected Rand Index e VI (Rank 1), mentre risulta meno performante rispetto ad altri metodi, ma compensato dagli altri punteggi, per il Dunn Index (Rank 21).

Visualizziamo i cluster ottenuti.

fviz_cluster(list(data = df, cluster = cluster_cut), geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

3.1.2. Dataset trasformato

Vediamo come cambia il clustering se invece che considerare tutte le variabili del dataset, consideriamo solo le tre componenti principali trovate prima.

banknote_pca_num <- banknote_pca[,-1]

results_all_pca <- evaluate_clustering(banknote_pca_num, dist(banknote_pca_num), status)

results_ranked_pca <- results_all_pca %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall_pca <- results_ranked_pca %>% arrange(total_rank) %>% head(1)
print(best_overall_pca)
##    Method k Silhouette      Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 average 2   0.431344 0.1056305      0.9020088 0.1997658        1        27
##   rank_corrected_rand rank_vi total_rank
## 1                   2       1         31

TODO: sintetizzare Il risultato indica che, applicando il clustering agglomerativo con numero di cluster pari a 2 e metodo di linkage “average” al dataset ridotto a tre componenti principali, otteniamo la combinazione migliore in base ai 4 indici di valutazione, anche se i risultati presentano alcune peculiarità. Rispetto alle altre combinazioni testate, la silhouette e la VI hanno ottenuto i migliori punteggi (rank 1), il Corrected Rand è quasi il migliore (rank 2), invece l’indice Dunn, in questa combinazione, non performa bene rispetto ad altre configurazioni testate (rank 27). Rispetto al clustering agglomerativo applicato su tutte le variabili del dataset il clustering sul dataset ridotto mostra: - Un indice silhouette leggermente migliore (0.431) rispetto a quello su tutte le variabili (0.375), suggerendo che i punti sono, mediamente, meglio raggruppati all’interno dei cluster nel caso PCA. - Il clustering su tutte le variabili ha un indice Dunn migliore (0.211 vs 0.106), il che significa che la separazione minima tra cluster è maggiore e i cluster risultano più distinti. - Entrambe le soluzioni mostrano valori molto alti del Corrected Rand e bassi del VI, ma la soluzione con tutte le variabili è leggermente superiore (0.960 vs 0.902 e 0.098 vs 0.200), indicando una migliore aderenza ad una struttura di riferimento o una maggiore stabilità del clustering.

Mentre l’approccio con PCA (riduzione a 3 componenti) offre cluster con una buona compattezza (indice silhouette migliore), il clustering basato su tutte le variabili, utilizzando il metodo ward.D2, risulta complessivamente superiore per quanto riguarda la separazione tra cluster e la corrispondenza con una struttura di riferimento (migliori valori di Dunn, Corrected Rand e VI).

Vediamo graficamente come sono separati i due cluster.

final_pca <- eclust(banknote_pca_num, "hclust", k = 2, hc_metric = "euclidean", hc_method = "average", graph = FALSE)

clusters <- final_pca$cluster

scatterplot3d(banknote_pca_num,
              color = clusters,
              pch = 19,
              main = "Clustering su PCA (PC1, PC2, PC3)",
              xlab = "PC1", ylab = "PC2", zlab = "PC3")

legend("topright", legend = unique(clusters), col = unique(clusters), pch = 19, title = "Cluster")

3.1.3. Dataset ridotto

Vediamo infine come cambia il clustering considerando solo le variabili Bottom e Diagonal che sembravano suddividere meglio i dati in gruppi.

banknote_bottom_diagonal <- banknote %>% select(Bottom, Diagonal)

banknote_bd_scaled <- scale(banknote_bottom_diagonal)

results_all_bd <- evaluate_clustering(banknote_bd_scaled, dist(banknote_bd_scaled), status)

results_ranked_bd <- results_all_bd %>%
  mutate(
    rank_sil = rank(-Silhouette, ties.method = "min"),
    rank_dunn = rank(-Dunn, ties.method = "min"),
    rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
    rank_vi = rank(VI, ties.method = "min"),
    
    total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
  )

best_overall_bd <- results_ranked_bd %>% arrange(total_rank) %>% head(1)
print(best_overall_bd)
##    Method k Silhouette       Dunn Corrected_Rand        VI rank_sil rank_dunn
## 1 ward.D2 2  0.6245096 0.09556396         0.9602 0.1120031        1        11
##   rank_corrected_rand rank_vi total_rank
## 1                   2       2         16

TODO: sintetizzare In questo caso la combinazione migliore è data da metodo di linkage ward.D2 e numero cluster pari a 2. Notiamo che in questo caso che il valore Silhouette 0.6245 è più alto rispetto agli altri due approcci (0.431 per PCA e 0.375 usando tutte le variabili), il che indica che i punti all’interno dei cluster sono ben raggruppati e distinti l’uno dall’altro. Ciò supporta l’idea iniziale che queste due variabili siano in grado di dividere bene il dataset. Il Dunn index risulta più basso rispetto a quello ottenuto con tutte le variabili (0.211) o con la PCA (0.106), il che potrebbe significare che alcuni cluster risultano più vicini o con una certa dispersione interna. Il Corrected Rand 0.9602 rimane molto elevato, segnalando una forte coerenza con la struttura di riferimento attesa. L’indice VI è basso, confermando una buona similarità o una discreta aderenza a una struttura di riferimento.

Visualizziamo come il dataset è stato diviso in cluster.

final_bd <- eclust(banknote_bd_scaled, "hclust", k = 2, hc_metric = "euclidean", hc_method = "ward.D2", graph = FALSE)

fviz_cluster(final_bd, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())

4. Clustering Partizionale

TODO: inserire commento introduttivo

4.1. K-Means

TODO: inserire commento introduttivo

4.1.1. Dataset completo

TODO: inserire commento introduttivo

4.1.1.1. Scelta del miglior numero di cluster k

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

df <- scale(banknote[, -1])

Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df, kmeans, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2)+
  labs(subtitle = "Elbow method")

Il metodo Silhouette.

fviz_nbclust(df, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

Il primo metodo e il secondo indicano 2 come numero ottimale di cluster, il terzo ne suggerisce 3. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

4.1.1.2. Validazione

TODO: inserire questa parte, prendere spunto dal clustering gerarchico

4.1.2 Dataset ridotto

4.1.2.1. Scelta del miglior numero di cluster k

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

df_ridotto <- df[, c("Bottom", "Diagonal")]

TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df_ridotto, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")

TODO: Il metodo Silhouette.

fviz_nbclust(df_ridotto, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

TODO: Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_ridotto, kmeans, nstart = 25, method = "gap_stat", nboot = 500)+
  labs(subtitle = "Gap statistic method")

TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

4.1.2.2. Validazione

TODO: da fare

4.1.3 Dataset trasformato

4.1.3.1. Scelta del miglior numero di cluster k

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.

df_pca <- scale(banknote_pca[, -1])

TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.

fviz_nbclust(df_pca, kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Elbow method")

TODO: Il metodo Silhouette.

fviz_nbclust(df_pca, kmeans, method = "silhouette")+
  labs(subtitle = "Silhouette method")

TODO: Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_pca, kmeans, nstart = 25, method = "gap_stat", nboot = 500, verbose=FALSE)+
  labs(subtitle = "Gap statistic method")

TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.

4.1.3.2. Validazione

TODO: da fare

4.2 K-Medoids

TODO: da rivedere tutto

4.2.1 Dataset completo

4.2.1.1. Scelta del miglior numero di cluster k

Come fatto in precedenza per l’algoritmo k-means, in questo capitolo si cercherà di determinare il numero di cluster k utilizzando i metodi elbow. silhouette e della GAP analysis sfruttando la funzione pam.

Vediamo il metodo Elbow.

fviz_nbclust(df, pam, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2)+
  labs(subtitle = "Elbow method")

Il metodo Silhouette.

fviz_nbclust(df, pam, method = "silhouette")+
  labs(subtitle = "Silhouette method")

Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df, pam, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

Tutti metodi suggeriscono che il numero ottimale dei cluster è pari a 2.

library(cluster)
head(df, n = 3)
##          Length      Left     Right     Bottom       Top  Diagonal
## [1,] -0.2549435  2.433346  2.829942 -0.2890067 -1.183765 0.4482473
## [2,] -0.7860757 -1.167507 -0.634788 -0.9120152 -1.432847 1.0557460
## [3,] -0.2549435 -1.167507 -0.634788 -0.4966762 -1.308306 1.4896737
pam.res <- pam(df, 2)
print(pam.res)
## Medoids:
##       ID     Length       Left      Right     Bottom        Top  Diagonal
## [1,] 193 -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,]  47 -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
## Clustering vector:
##   [1] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Objective function:
##    build     swap 
## 1.842553 1.788867 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       2
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##          Length       Left      Right     Bottom        Top  Diagonal
## [1,] -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,] -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
head(pam.res$clustering)
## [1] 1 2 2 2 2 2
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

4.2.2 Dataset ridotto

4.2.2.1. Scelta del miglior numero di cluster k

Come fatto in precedenza per l’algoritmo k-means, in questo capitolo si cercherà di determinare il numero di cluster k utilizzando i metodi elbow. silhouette e della GAP analysis sfruttando la funzione pam.

Vediamo il metodo Elbow.

df_ridotto <- df[, c("Bottom", "Diagonal")]
fviz_nbclust(df_ridotto, pam, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2)+
  labs(subtitle = "Elbow method")

Il metodo Silhouette.

fviz_nbclust(df_ridotto, pam, method = "silhouette")+
  labs(subtitle = "Silhouette method")

Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_ridotto, pam, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

df_ridotto <- df[, c("Bottom", "Diagonal")]
pam.res <- pam(df_ridotto, 2)
print(pam.res)
## Medoids:
##       ID     Bottom   Diagonal
## [1,]  47 -0.7735689  0.8821750
## [2,] 185  0.8877871 -0.8535357
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Objective function:
##     build      swap 
## 0.8831314 0.6414212 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       1
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       1
pam.res$medoids
##          Bottom   Diagonal
## [1,] -0.7735689  0.8821750
## [2,]  0.8877871 -0.8535357
head(pam.res$clustering)
## [1] 1 1 1 1 1 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

4.2.3 Dataset trasformato

4.2.3.1. Scelta del miglior numero di cluster k

Come fatto in precedenza per l’algoritmo k-means, in questo capitolo si cercherà di determinare il numero di cluster k utilizzando i metodi elbow. silhouette e della GAP analysis sfruttando la funzione pam.

Vediamo il metodo Elbow.

df_pca <- scale(banknote_pca[, -1])
fviz_nbclust(df_ridotto, pam, method = "wss") +
  geom_vline(xintercept = 2, linetype = 2)+
  labs(subtitle = "Elbow method")

Il metodo Silhouette.

fviz_nbclust(df_pca, pam, method = "silhouette")+
  labs(subtitle = "Silhouette method")

Ed infine la GAP analysis:

set.seed(123)
fviz_nbclust(df_pca, pam, nstart = 25, method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

pam.res <- pam(df_pca, 2)
print(pam.res)
## Medoids:
##      ID       PCA1       PCA2       PCA3
## 198 198  0.8816755 -0.1289878 0.19323058
## 21   21 -1.1983358  0.0749874 0.05222386
## Clustering vector:
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   2   2   2   2   1   2   2   2   1   1   2   2   2   2   2   2   2   2   2 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   2   2   2   2   2   1   2   2   2   1   2   2   2   2   2   2   2   2   2   2 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## Objective function:
##    build     swap 
## 1.370309 1.329691 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##           PCA1       PCA2       PCA3
## 198  0.8816755 -0.1289878 0.19323058
## 21  -1.1983358  0.0749874 0.05222386
head(pam.res$clustering)
## 1 2 3 4 5 6 
## 1 2 2 2 2 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

TODO: da rivedere tutto

4.2.1 Dataset completo

In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-medoids con k impostato a 2.

Come in precedenza, si analizzano prima i grafici a dispersione e successivamente i due cluster.

head(df, n = 3)
##          Length      Left     Right     Bottom       Top  Diagonal
## [1,] -0.2549435  2.433346  2.829942 -0.2890067 -1.183765 0.4482473
## [2,] -0.7860757 -1.167507 -0.634788 -0.9120152 -1.432847 1.0557460
## [3,] -0.2549435 -1.167507 -0.634788 -0.4966762 -1.308306 1.4896737
pam.res <- pam(df, 2)
print(pam.res)
## Medoids:
##       ID     Length       Left      Right     Bottom        Top  Diagonal
## [1,] 193 -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,]  47 -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
## Clustering vector:
##   [1] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Objective function:
##    build     swap 
## 1.842553 1.788867 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       2
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##          Length       Left      Right     Bottom        Top  Diagonal
## [1,] -0.5205096  0.4944248  0.6026155  0.9570103  0.5598130 -1.113892
## [2,] -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584  0.882175
head(pam.res$clustering)
## [1] 1 2 2 2 2 2
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

4.2.2 Dataset ridotto

df_ridotto <- df[, c("Bottom", "Diagonal")]
pam.res <- pam(df_ridotto, 2)
print(pam.res)
## Medoids:
##       ID     Bottom   Diagonal
## [1,]  47 -0.7735689  0.8821750
## [2,] 185  0.8877871 -0.8535357
## Clustering vector:
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Objective function:
##     build      swap 
## 0.8831314 0.6414212 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       1
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       1
pam.res$medoids
##          Bottom   Diagonal
## [1,] -0.7735689  0.8821750
## [2,]  0.8877871 -0.8535357
head(pam.res$clustering)
## [1] 1 1 1 1 1 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

4.2.3 Dataset trasformato

df_pca <- scale(banknote_pca[, -1])
pam.res <- pam(df_pca, 2)
print(pam.res)
## Medoids:
##      ID       PCA1       PCA2       PCA3
## 198 198  0.8816755 -0.1289878 0.19323058
## 21   21 -1.1983358  0.0749874 0.05222386
## Clustering vector:
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   2   2   2   2   1   2   2   2   1   1   2   2   2   2   2   2   2   2   2 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   2   2   2   2   2   1   2   2   2   1   2   2   2   2   2   2   2   2   2   2 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## Objective function:
##    build     swap 
## 1.370309 1.329691 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)
##    Status Length  Left Right Bottom  Top Diagonal cluster
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0       1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7       2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2       2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0       2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8       2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4       1
## 7 genuine  215.5 129.5 129.7    7.9  9.6    141.6       2
## 8 genuine  214.5 129.6 129.2    7.2 10.7    141.7       2
pam.res$medoids
##           PCA1       PCA2       PCA3
## 198  0.8816755 -0.1289878 0.19323058
## 21  -1.1983358  0.0749874 0.05222386
head(pam.res$clustering)
## 1 2 3 4 5 6 
## 1 2 2 2 2 1
cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])

Come previsto, i cluster risultano molto simili.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic()
)

5. Clustering Model-based: Mixture of Gaussian

In questa sezione introduciamo l’applicazione del clustering basato su misture di gaussiane, utilizzando la funzione Mclust del pacchetto mclust. Questo approccio assume che i dati siano generati da una combinazione di distribuzioni gaussiane e, tramite l’algoritmo EM (Expectation-Maximization), stima i parametri di ciascuna componente. Il vantaggio principale di questo metodo è la capacità di modellare cluster con forme ellissoidali, permettendo di gestire anche cluster con sovrapposizioni. Inoltre, Mclust effettua una selezione automatica del numero ottimale di cluster e del modello di covarianza (parsimonioso) in base al criterio BIC, fornendo un’indicazione della bontà dell’adattamento e della complessità del modello.

data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 4 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -1191.595 200 51 -2653.405 -2666.898
## 
## Clustering table:
##  1  2  3  4 
## 16 99 47 38
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

Nonostante sappiamo che il dataset banknote presenta due classi reali (ad esempio, banconote fraudolente e non fraudolente), l’analisi con Mclust ha individuato come miglior modello un numero di cluster pari a 3. Inoltre, il modello selezionato è di tipo VEE (ovvero ellissoidale, con forma variabile, ma con volume ed orientamento uguali), in base al valore più alto di BIC evidenziato nel plot. È interessante notare come il BIC per 2 cluster risulti significativamente inferiore rispetto a quello per un numero maggiore di cluster, probabilmente perché sono state utilizzate tutte le variabili del dataset per costruire i cluster.
Confrontiamo esplicitamente i valori di BIC e ICL per soluzioni con 2 e 3 cluster:

data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -1260.823 200 49 -2781.264 -2781.417
## 
## Clustering table:
##   1   2 
## 101  99
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: 
## 
##  log-likelihood   n df      BIC      ICL
##       -1212.796 200 43 -2653.42 -2653.74
## 
## Clustering table:
##  1  2  3 
## 16 99 85

Dai sommari ottenuti si osserva che entrambi i modelli, per k = 2 e k = 3, forniscono valori che suggeriscono una soluzione migliore per 3 cluster. I valori del BIC e, conseguentemente, anche quelli dell’ICL (che include una penalizzazione per l’incertezza nelle assegnazioni), indicano che la soluzione a 3 cluster è preferibile rispetto a quella a 2 cluster.
Valutiamo l’effetto della scelta delle variabili sul modello. Utilizziamo invece solo le due variabili “Bottom” e “Diagonal”, individuate in precedenza dagli scatterplot come le più discriminanti:

data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -387.8563 200 13 -844.5908 -849.2887
## 
## Clustering table:
##   1   2   3 
## 100  16  84
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

Anche in questo caso si nota che il valore di BIC è notevolmente più alto (migliore) rispetto al caso in cui sono state usate tutte le variabili. L’algoritmo suggerisce ancora un numero ottimale di cluster pari a 3, ma questa volta ha selezionato come modello parsimonioso il modello EVE (ellissoidale, con volume uguale e orientamento variabile). È importante osservare che in uno dei cluster sono presenti solamente 16 osservazioni, il che potrebbe indicare un gruppo di dati che, probabilmente, sono state erroneamente assegnate a un cluster separato.
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC confrontando esplicitamente le soluzioni per 2 e 3 cluster:

data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVE (ellipsoidal, equal orientation) model with 2 components: 
## 
##  log-likelihood   n df      BIC       ICL
##       -417.7539 200 10 -888.491 -890.1796
## 
## Clustering table:
##   1   2 
##  99 101
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -387.8563 200 13 -844.5908 -849.2887
## 
## Clustering table:
##   1   2   3 
## 100  16  84

In questo caso, i valori di BIC e ICL sono concordi: per un numero di cluster pari a 3 si ottiene un BIC di circa -845 e un ICL di -849, mentre per 2 cluster il BIC risulta attorno a -888 e l’ICL a -890. Questo dimostra che anche l’ICL, che penalizza ulteriormente le assegnazioni incerte, suggerisce che il modello ottimale è quello con 3 cluster.
Adesso, procederemo ad utilizzare solamente le variabili estratte tramite la tecnica PCA.

dati <- banknote_pca[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -764.5395 200 16 -1613.852 -1616.362
## 
## Clustering table:
##   1   2 
## 102  98
plot(gmm_result, what = "BIC")

plot(gmm_result, what = "classification")

plot(gmm_result, what = "uncertainty")

plot(gmm_result, what = "density")

In questo caso si nota che il valore di BIC è più alto (migliore) rispetto al caso in cui sono state usate tutte le variabili ma è più basso del caso in cui sono state usate solamente le due variabili selezionate. Questa volta l’algoritmo suggerisce un numero ottimale di cluster pari a 2, selezionando come modello parsimonioso il modello EEV (ellissoidale, con volume e orientamento uguale).
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC confrontando esplicitamente le soluzioni per 2 e 3 cluster:

gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -764.5395 200 16 -1613.852 -1616.362
## 
## Clustering table:
##   1   2 
## 102  98
gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##        -764.778 200 19 -1630.224 -1657.822
## 
## Clustering table:
##  1  2  3 
## 64 98 38

Anche in questo caso, i valori di BIC e ICL sono concordi: per un numero di cluster pari a 3 si ottiene un BIC di circa -1630 e un ICL di -1657, mentre per 2 cluster il BIC risulta attorno a -1614 e l’ICL a -1616. Questo dimostra che anche l’ICL suggerisce che il modello ottimale è quello con 2 cluster.

6. Classificazione con gli algoritmi di cluster

In questa sezione per le tre tipologie di dataset (completo, ridotto e trasformato), analizzeremo la qualità dei diversi metodi di clustering, fissando il numero di cluster pari al numero di classi reali, ovvero 2. Questo approccio ci permette di utilizzare tecniche di validazione esterne, dato che disponiamo delle etichette reali delle osservazioni (variabile Status).
Per valutare l’efficacia dei cluster identificati, utilizzeremo diverse metriche, tra cui:
- Adjusted Rand Index (ARI) che misura il grado di concordanza tra il clustering ottenuto e la classificazione reale. - Metriche di valutazione proprie del machine learning supervisionato, come accuratezza, specificità, recall, matrice di confusione e AUC, per interpretare il clustering come un problema di classificazione binaria.

6.1. Dataset completo

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset originale con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo. Da precisare che in aggiunta agli algoritmi discussi in precedenza, confrontiamo anche l’algoritmo gerarchico divisivo che non è stato affrontato prima.

data("banknote")
data <- banknote[, -1]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      1        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      2        2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      2        2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      2        2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      2        2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      1        2
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   2
## 2                          2                     2   2
## 3                          2                     2   2
## 4                          2                     2   2
## 5                          2                     2   2
## 6                          2                     2   2

Il codice successivo implementa una funzione per riallineare le assegnazioni dei cluster alle classi reali in modo da rendere i confronti più accurati.

bestMapping <- function(true, pred) {

  banknote <- banknote %>%
  mutate(Status_numeric = recode(Status,
                                 "counterfeit" = 1,
                                 "genuine" = 2))


  true <- banknote$Status_numeric
  true <- as.character(true)
  pred <- as.character(pred)
  levels_true <- sort(unique(true))
  
  cm1 <- confusionMatrix(factor(pred, levels = levels_true), factor(true, levels = levels_true))
  acc1 <- cm1$overall["Accuracy"]
  
  pred_swap <- ifelse(pred == levels_true[1], levels_true[2], levels_true[1])
  cm2 <- confusionMatrix(factor(pred_swap, levels = levels_true), factor(true, levels = levels_true))
  acc2 <- cm2$overall["Accuracy"]
  
  if(acc2 > acc1) {
    return(list(mapped = pred_swap, cm = cm2))
  } else {
    return(list(mapped = pred, cm = cm1))
  }
}

Nel blocco successivo, calcoliamo le metriche di valutazione per ciascun algoritmo di clustering. È importante notare che il valore di AUC ha un significato più chiaro per l’algoritmo GMM (Gaussian Mixture Model), poiché questo metodo fornisce direttamente le probabilità a posteriori di appartenenza ai cluster. Per gli altri algoritmi, che operano in modo hard (assegnando ciascun punto a un solo cluster senza una misura di incertezza), il calcolo dell’AUC è stato forzato attraverso:

  • Per K-means e K-medoids, è stato definito uno score basato sulla distanza dei punti dal centroide o dal medoide del rispettivo cluster.
  • Per gli algoritmi gerarchici, non è stato possibile adottare un criterio analogo per costruire uno score continuo, motivo per cui l’AUC non è stato calcolato (indicata come NA nella tabella).
# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.8456292    0.960        1.00        0.92 0.9973
## 2                K-medoids 0.9406018    0.985        1.00        0.97 0.9991
## 3 Gerarchico Agglomerativo 0.9602001    0.990        1.00        0.98     NA
## 4      Gerarchico Divisivo 0.9406015    0.985        0.99        0.98     NA
## 5                      GMM 0.9799995    0.995        1.00        0.99 0.9999

I risultati mostrano che il metodo GMM ha la migliore performance complessiva, con il più alto valore di ARI (0.98) e un’accuratezza del 99.5%. Il clustering gerarchico agglomerativo ha anch’esso buone prestazioni, mentre K-means e K-medoids mostrano una buona affidabilità ma con AUC leggermente inferiori. Analizzando la sensibilità e la specificità, si nota che la sensibilità – ovvero la capacità di identificare correttamente le banconote contraffatte – è generalmente molto alta per tutti gli algoritmi, indicando un’elevata capacità di individuare i falsi. Tuttavia, la specificità, ossia la capacità di riconoscere correttamente le banconote genuine, varia maggiormente tra i metodi. In particolare, K-means mostra una specificità inferiore (0.92) rispetto a K-medoids (0.97) e ai metodi gerarchici (0.98), suggerendo una maggiore tendenza a classificare erroneamente alcune banconote genuine come contraffatte; questo errore è meno grave di classificare una banconota contraffatta come vera.

6.2. Dataset ridotto

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset ridotto (con le sole variabili Diagonal e Bottom) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.

data("banknote")
data <- banknote[, c("Bottom", "Diagonal")]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      2        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      2        1
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      2        1
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      2        1
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      2        1
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      2        1
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   1
## 2                          1                     1   1
## 3                          1                     1   1
## 4                          1                     1   1
## 5                          1                     1   1
## 6                          1                     1   1

Nel blocco successivo, calcoliamo le metriche di valutazione come nel caso precedente.

# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.9799995    0.995        0.99        1.00 0.9999
## 2                K-medoids 0.9799995    0.995        0.99        1.00 0.9999
## 3 Gerarchico Agglomerativo 0.9602000    0.990        0.99        0.99     NA
## 4      Gerarchico Divisivo 0.9406015    0.985        0.98        0.99     NA
## 5                      GMM 0.9799995    0.995        1.00        0.99 0.9995

Questa configurazione semplifica il modello senza perdere accuratezza, dimostrando che Bottom e Diagonal sono le variabili chiave per distinguere le classi:
- K-means, K-medoids e GMM ottengono ARI 0.98, accuracy 99.5% e AUC ~1, confermando che queste due variabili sono altamente discriminative.
- GMM raggiunge sensibilità 100%, identificando tutte le banconote contraffatte, mentre K-means e K-medoids garantiscono specificità 100%, evitando falsi positivi.

6.3. Dataset ridotto

Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset trasformato (con le variabili trasformate mediante tecnica PCA) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.

data <- banknote_pca[, -1]
data_scaled <- scale(data)

set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)

banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification

head(banknote)
##    Status Length  Left Right Bottom  Top Diagonal kmeans kmedoids
## 1 genuine  214.8 131.0 131.1    9.0  9.7    141.0      2        1
## 2 genuine  214.6 129.7 129.7    8.1  9.5    141.7      1        2
## 3 genuine  214.8 129.7 129.7    8.7  9.6    142.2      1        2
## 4 genuine  214.8 129.7 129.6    7.5 10.4    142.0      1        2
## 5 genuine  215.0 129.6 129.7   10.4  7.7    141.8      1        2
## 6 genuine  215.7 130.8 130.5    9.0 10.1    141.4      1        1
##   hierarchical_agglomerative hierarchical_divisive gmm
## 1                          1                     1   1
## 2                          2                     2   2
## 3                          2                     2   2
## 4                          2                     2   2
## 5                          2                     2   2
## 6                          1                     2   2

Nel blocco successivo, calcoliamo le metriche di valutazione come nei casi precedenti.

# K-means
ari_kmeans   <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans   <- cm_kmeans$overall["Accuracy"]
sens_kmeans  <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans  <- cm_kmeans$byClass["Specificity"]

distances_km <- t(apply(data_scaled, 1, function(x) {
  apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)

# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids   <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids  <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids  <- cm_kmedoids$byClass["Specificity"]

medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
  apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)

# Gerarchico Agglomerativo
ari_hierAgg  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg   <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg  <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg  <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA

# Gerarchico Divisivo
ari_hierDiv  <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv   <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv  <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv  <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA

# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm   <- cm_gmm$overall["Accuracy"]
sens_gmm  <- cm_gmm$byClass["Sensitivity"]
spec_gmm  <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)

Infine, raccogliamo i risultati in un dataframe riassuntivo.

results_df <- data.frame(
  Metodo       = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
  ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
  Accuracy     = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
  Sensitivity  = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
  Specificity  = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
  AUC          = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)
##                     Metodo       ARI Accuracy Sensitivity Specificity    AUC
## 1                  K-means 0.8830121    0.970        0.97        0.97 0.9987
## 2                K-medoids 0.8272389    0.955        1.00        0.91 0.9946
## 3 Gerarchico Agglomerativo 0.5160469    0.860        0.95        0.77     NA
## 4      Gerarchico Divisivo 0.9212040    0.980        0.98        0.98     NA
## 5                      GMM 0.9212042    0.980        0.99        0.97 0.9983

Rispetto ai risultati con il dataset completo o ridotto, l’uso di solo 3 variabili ottenute con PCA ha ridotto la performance globale, in particolare per K-medoids e Gerarchico Agglomerativo. Sebbene l’accuratezza resti elevata, l’ARI e la specificità mostrano una leggera diminuzione, il che suggerisce che la riduzione delle variabili ha compromesso la capacità di questi algoritmi di catturare la struttura dei dati. In particolare, l’algoritmo Gerarchico Agglomerativo mostra una performance molto più bassa in termini di ARI e accuratezza, indicando che l’algoritmo è più sensibile alla perdita di variabili discriminanti.

In conclusione, dal punto di vista della misurazione delle performance con tecniche esterne (avendo a disposizione la vera classe delle osservazioni), l’algoritmo GMM è risultato il migliore nei tre diversi dataset, anche se, tutto sommato, gli altri algoritmi non hanno ottenuto dei cattivi risultati.

7. Conclusioni

Questo progetto, realizzato a scopo didattico, ha analizzato diverse tecniche di clustering applicate al dataset banknote, un dataset semplice con due classi ben distinte. Inizialmente sono state utilizzate misure interne – quali il BIC, l’indice silhouette e altre metriche di compattezza – per valutare la qualità dei cluster senza ricorrere alle etichette reali. Successivamente, l’analisi è stata estesa con misure esterne (Adjusted Rand Index, accuratezza, sensibilità, specificità e AUC) per confrontare le assegnazioni dei cluster con la vera classe delle osservazioni.

I risultati mostrano che, quando si utilizzano tutte le variabili o solo quelle discriminanti (Bottom e Diagonal), algoritmi come GMM e, in alcuni casi, K-means e K-medoids raggiungono prestazioni eccellenti, con ARI e AUC quasi perfetti. Al contrario, i metodi gerarchici, in particolare quello agglomerativo, hanno ottenuto risultati leggermente inferiori, anche se in alcune configurazioni il clustering gerarchico divisivo si è dimostrato competitivo. L’approccio GMM risulta particolarmente robusto, grazie alla possibilità di estrarre score continui tramite le probabilità a posteriori.

Questi risultati, sia interni che esterni, evidenziano come una corretta selezione o trasformazione delle variabili (ad es. tramite PCA) possa semplificare il modello senza compromettere la capacità di distinguere correttamente le banconote genuine da quelle contraffatte.

Per il futuro, un possibile sviluppo potrebbe essere l’applicazione di questi metodi a dataset più complessi, nonché l’ottimizzazione degli score per metodi hard. Inoltre, integrare tecniche di soft clustering o combinare il clustering con approcci supervisionati potrebbe migliorare ulteriormente la generalizzazione e la robustezza dei modelli, offrendo una visione più completa delle performance in contesti reali.